| Weights matrix used | Global Moran’s I value |
|---|---|
| contiguity, 1st order | 0.472 |
| contiguity, 2nd order | 0.212 |
| distance-based | 0.283 |
All of these Moran’s I values are postive, indicating that the variable shows positive spatial autocorrelation. The Moran’s I value changes depending on the weights matrix used. The Moran’s I value calculated using the 1st order contiguity weights matrix is the most postive. This weights matrix considers all contiguous polygons using queens adjacency to be neighbors. The 2nd order contiguity takes neighbors of neighbors, so it makes sense that using this weights matrix we see less strong positive autocorrelation. The distance based weights matrix takes the distance between centroids to determine adjacency – this ends up giving fewer neighbors to larger polygons. This may have something to do with its relatively small Moran’s I value compared to the contiguity, 1st order value.
| Weights matrix used | # significant Moran’s I | # in high-high clusters | # in low-low clusters |
|---|---|---|---|
| contiguity, 1st order | 7 | 16 | 8 |
| contiguity, 2nd order | 9 | 17 | 10 |
| distance-based | 8 | 11 | 13 |
| Weights matrix used | # significant Gi | # in high clusters | # in low clusters |
|---|---|---|---|
| contiguity, 1st order | 7 | 17 | 9 |
| contiguity, 2nd order | 7 | 19 | 16 |
| distance-based | 7 | 11 | 17 |
library(rgdal)
## Loading required package: sp
## rgdal: version: 1.2-5, (SVN revision 648)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.1.2, released 2016/10/24
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
## Linking to sp version: 1.2-4
datapath = "/users/miadawson/Documents/~UC Davis~/~GEO 200CN/labs/lab 9"
nepal=readOGR(dsn=datapath, layer="Nepal")
## OGR data source with driver: ESRI Shapefile
## Source: "/users/miadawson/Documents/~UC Davis~/~GEO 200CN/labs/lab 9", layer: "Nepal"
## with 75 features
## It has 26 fields
## Integer64 fields read as strings: ID KIDS1_5 AD_ILGT50
hist(nepal$MALKIDS, main='Rate of childhood malnutrition in Nepalese districts', xlab = NULL)
library (tmap)
#vignette(package = "tmap") # available vignettes in tmap
#vignette("tmap-nutshell") #gives you an overview of the package
tmap_mode('plot')
## tmap mode set to plotting
qtm (shp=nepal, fill="MALKIDS", fill.palette="Blues") #map just the MALKIDS
qtm (shp=nepal, fill=c("MALKIDS", "POPULATION"), fill.palette="Blues") #map both
tmap has two modes “plot” return static maps and “view” returns interactive maps. You toggle between them with the function tmap_mode() [ tmap_mode(”view”) or tmap_mode(”plot”)]. Switch to “view” mode and plot the maps above again.
tmap_mode('view')
## tmap mode set to interactive viewing
qtm (shp=nepal, fill="MALKIDS", fill.palette="Blues") #map just the MALKIDS
qtm (shp=nepal, fill=c("MALKIDS", "POPULATION"), fill.palette="Blues") #map both
Now let’s plot a scatterplot similar to what we saw in GeoDa.
library (ggplot2)
library(labeling)
ggplot(nepal@data, aes(POPULATION, MALKIDS)) + geom_point(aes(col=MALKIDS, size=MALKIDS))
Both tmap and ggplot2 work by adding layers so in this case we set the general plot tell it that the data comes from npl@data, and the values for the aesthetic mapping come from population, and malkids. Then we add another layer, the actual points where the color and size of the points are defined by the value of malkids.
Looking at the polt you might decide you want the higher values to have the darker colors, you can do that by simply adding a negative sign in front of MALKIDS [col=–MALKIDS].
library(raster)
library(spdep)
## Loading required package: Matrix
##1st order queen cont
?poly2nb
w <- poly2nb(nepal, row.names=nepal$ID) ##default is queen=TRUE
ww <- nb2listw(w, style='B')
##2nd order queen cont
w2 = nblag(w, 2)
w22=w2[[2]] #subset for only 2nd order neighbors
ww22 <- nb2listw(w22, style='B')
##distance-based cont
dw = poly2nb(nepal, row.names=nepal$ID, snap = 0.7)
dww <- nb2listw(dw, style='B')
##1st order queen cont
moran(nepal$MALKIDS, ww, n=length(ww$neighbours), S0=Szero(ww))
## $I
## [1] 0.4798273
##
## $K
## [1] 2.430716
##2nd order queen cont
moran(nepal$MALKIDS, ww22, n=length(ww22$neighbours), S0=Szero(ww22))
## $I
## [1] 0.2258045
##
## $K
## [1] 2.430716
##distance-based cont
moran(nepal$MALKIDS, dww, n=length(dww$neighbours), S0=Szero(dww))
## $I
## [1] 0.2616114
##
## $K
## [1] 2.430716
| Weights matrix used | Moran’s I value (GeoDa) | Moran’s I Value (R) |
|---|---|---|
| contiguity, 1st order | 0.472 | 0.479 |
| contiguity, 2nd order | 0.212 | 0.225 |
| distance-based | 0.283 | 0.261 |
These values are close, but not exactly the same as the values I got in GeoDa. There may be slight differences in the way GeoDa and the R package calculate these values. For example, if the programs round differently within a complicated expression, this could lead to different output values for Moran’s I.
Now let’s move into looking at some of the local patterns. First we’ll calculate the Local Moran’s I, and create a dataframe with the values, lag values and p-vlaues to then we visulaize the results in similar ways to what we saw in GeoDa.
Again test with each of your spatial weight matrcies.
##1st order queens
LoM <- localmoran(nepal$MALKIDS, ww, zero.policy=NULL)
LoMdf <- data.frame (LoM)
lag = lag.listw(ww, nepal$MALKIDS)
dat = data.frame(ID = nepal$ID, mal = nepal$MALKIDS, lag = lag, pval = LoMdf$Pr.z...0)
dat$sig =ifelse (dat$pval < 0.05 ,1, 0 ) #set signifcance at p=0.05
mv = mean(dat$mal) #create a mean for the value
mh = mean (dat$lag) #create a mean for the lagged value
ggplot(dat, aes(x = mal, y = lag, col=-sig)) + geom_point() + geom_smooth(method = "lm", col="grey") +
geom_vline(xintercept=mv, linetype = "dashed", col="grey") + geom_hline(yintercept = mh, linetype="dashed",
col="grey")
##2nd order queens
LoM <- localmoran(nepal$MALKIDS, ww22, zero.policy=NULL)
LoMdf <- data.frame (LoM)
lag = lag.listw(ww22, nepal$MALKIDS)
dat = data.frame(ID = nepal$ID, mal = nepal$MALKIDS, lag = lag, pval = LoMdf$Pr.z...0)
dat$sig =ifelse (dat$pval < 0.05 ,1, 0 ) #set signifcance at p=0.05
mv = mean(dat$mal) #create a mean for the value
mh = mean (dat$lag) #create a mean for the lagged value
ggplot(dat, aes(x = mal, y = lag, col=-sig)) + geom_point() + geom_smooth(method = "lm", col="grey") +
geom_vline(xintercept=mv, linetype = "dashed", col="grey") + geom_hline(yintercept = mh, linetype="dashed",
col="grey")
##distance-based
LoM <- localmoran(nepal$MALKIDS, dww, zero.policy=NULL)
LoMdf <- data.frame (LoM)
lag = lag.listw(dww, nepal$MALKIDS)
dat = data.frame(ID = nepal$ID, mal = nepal$MALKIDS, lag = lag, pval = LoMdf$Pr.z...0)
dat$sig =ifelse (dat$pval < 0.05 ,1, 0 ) #set signifcance at p=0.05
mv = mean(dat$mal) #create a mean for the value
mh = mean (dat$lag) #create a mean for the lagged value
ggplot(dat, aes(x = mal, y = lag, col=-sig)) + #sets x axis as malnutrition rate, y value as lagged value, color as significance
geom_point() + #specifies scatterplot
geom_smooth(method = "lm", col="grey") + #adds smoothed line of best fit
geom_vline(xintercept=mv, linetype = "dashed", col="grey") + #adds vertical line bisecting data
geom_hline(yintercept = mh, linetype="dashed", col="grey") #adds horizontal line bisecting data
Try to unpick this rather long line, and see if you understand the different layers that are making up this plot. Looking at this plot what can you say about the pattern of Local Moran’s I? (see above comments)
From these plots, it appears that there is most postive spatial autocorrelation when using the 1st order queens contiguity. According to OSU, positive spatial autocorrelation is evidenced by more data falling in the bottom left and upper right corners of Moran’s scatterplots. This pattern is most visibile in the plot using 1st order quuen’s contiguity. So, it should have the highest local Moran’s I. The 2nd order queen’s contiguity based plot appears to have a less positive local moran’s I, and the distance based seems to have the least positive local moran’s I.
For the Gi and Gi* go read and follow along with the with the Rspatial chapter on local Statistics to create local statistics for the Nepal dataset[http://rspatial.org/rosu/rst/Chapter8.html], stop with the Local Moran’s I ie don’t work through the Geographically Weighted Regression section.
##local getis Gi
wr <- poly2nb(nepal, row.names=nepal$ID, queen=FALSE)
lstw <- nb2listw(wr, style='B')
Gi <- localG(nepal$MALKIDS, lstw)
head(Gi)
## [1] 0.4309344 1.5261740 0.2739212 -0.8286501 -3.2037288 -1.7890001
#plot
par(mai=c(0,0,0,0))
Gcuts <- cut(Gi, 5)
Gcutsi <- as.integer(Gcuts)
cols <- rev(gray(seq(0,1,.2)))
plot(nepal, col=cols[Gcutsi])
legend('bottomleft', levels(Gcuts), fill=cols)
##classify each polygon as own neighbor
ws <- include.self(wr)
lstws <- nb2listw(ws, style='B')
Gis <- localG(nepal$MALKIDS, lstws)
Gscuts <- cut(Gis, 5)
Gscutsi <- as.integer(Gscuts)
cols <- rev(gray(seq(0,1,.2)))
plot(nepal, col=cols[Gscutsi])
legend('bottomleft', levels(Gscuts), fill=cols)
#local average
m <- sapply(ws, function(i) mean(nepal$MALKIDS[i]))
cts <- cut(m, 5)
mcts <- as.integer(cts)
plot(nepal, col=cols[mcts])
legend('bottomleft', levels(cts), fill=cols)
#local Moran's i
Ii <- localmoran(nepal$MALKIDS, lstw)
Icuts <- cut(Ii, 5)
Icutsi <- as.integer(Icuts)
plot(nepal, col=cols[Icutsi])
legend('bottomleft', levels(Icuts), fill=cols)
These values are similar to what I got in GeoDa. As before, there may be differences in the way GeoDa and the R package calculate these values. For example, if the programs round differently within a complicated expression, this could lead to different output values.